function [VaR] = VaRCov(rPortfolio,rIndex,beta,weights, alpha, V0)

% Note that rPortfolio is a vector of return (rPortfolio = rPortfolioortfolio * weights).
% rPortfolioortfolio is a matrix of return of stocks.
% ----------------------------------------------------    

    % VaR of Market
    varMarket = var(rIndex);
      
    % Create a vector of error variance
    [lenStock,numStock] = size(rPortfolio);
    %vec_evar = cov(rPortfolio(:,1),1);
    rPortfolio_mean = mean(rPortfolio,1);
    rPortfolio_error = zeros(size(rPortfolio));
    for i=1:numStock
        rPortfolio_error(:,i) = rPortfolio(:,i) - rPortfolio_mean(i);
        D(i) = var(rPortfolio_error(:,i));
    end
    
    D = diag(D);
    
%     vec_evar = cov(rPortfolio_error);%[vec_evar      cov(rPortfolio(:,i),1)];
%    
%     
%     % Diagonal matrix of error variance
%     D = diag(vec_evar(1:size(vec_evar,1)+1:end));

    

    % z-statistic that corresponds to the 1-alpha percentile of a standard normal distribution.
    Z = norminv(alpha, 0, 1); 
    
    %% Calculate the covariance Matrix Sigma
    % Consider 2 cases: beta is a constant and beta is a vector
    
    m = size(beta);
    beta = beta';
    
    if m(1) == 1
        SigmaT = beta * beta' * varMarket + D;
        
        % Calculate the VaR of Portfolio 
        VaR = -Z * sqrt(weights' * SigmaT * weights) * V0;  
    else
       for i=1:m(1)
           SigmaT = beta(:,i) * beta(:,i)' * varMarket + D;
           VaR(i) = -Z * sqrt(weights' * SigmaT * weights) * V0;  
       end
    end
    VaR = VaR';